1 Set up

Install necessary packages:

install.packages("tidyverse", dependencies = TRUE) 
install.packages("raster", dependencies = TRUE) 
install.packages("sf", dependencies = TRUE) 
install.packages("landscapetools", dependencies = TRUE)
install.packages("landscapemetrics", dependencies = TRUE)

2 Prepare data

2.1 Quick visualisation

Import only the RGB color bands as individual RasterLayer objects:

library(raster)

#blue
b2 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B2.tif')

#green
b3 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B3.tif')

#red
b4 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B4.tif')

Combine the RasterLayer objects and visualise the satellite image:

landsatRGB <- stack(b4, b3, b2) #order is impt

plotRGB(landsatRGB, 
        stretch = "lin") #scale the values (try using "hist" also)

Landsat-8 true color composite (RGB). Source: U.S. Geological Survey.

2.2 Import data

Import all 5 bands from the satellite data as a RasterStack object named landsat:

filenames <- paste0('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B', 1:5, ".tif")

landsat <- stack(filenames)

#rename bands
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR')

Check coordinate reference system of landsat:

crs(landsat)
## CRS arguments:
##  +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

2.3 Crop data

Import polygon of city boundaries as sgshp and check if their coordinate reference systems match:

library(sf)

sgshp <- st_read("data/master-plan-2014-region-boundary-web-shp/MP14_REGION_WEB_PL.shp")

Check coordinate reference system of sgshp:

crs(sgshp)
## [1] "+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs"

Transform sgshp to the match the coordinate reference system of the landsat:

sgshp <- st_transform(sgshp, crs = crs(landsat))

Crop landsat according to city boundaries:

landsat <- crop(landsat, sgshp) #crop to rectangle

landsat <- mask(landsat, sgshp) #mask values according to shape of sgshp

Plot the cropped image using only the RGB bands:

landsatRGB <- subset(landsat, c(4,3,2)) #Red, Green, Blue

plotRGB(landsatRGB,
        stretch = "lin")

Landsat-8 true color composite (USGS, 2018) cropped to city boundaries (URA, 2014)

3 Classify land cover

3.1 Calculate NDVI

Create a function that calcuates the Normalized Difference Vegetation Index (NDVI) for each pixel:

ndvi <- function(x, y) {
    (x - y) / (x + y)
  }

Apply function to the NIR and Red bands of landsat

landsatNDVI <- overlay(landsat[[5]], landsat[[4]], 
                       fun = ndvi)

Limit the range of values to be from -1 to 1:

landsatNDVI <- reclassify(landsatNDVI, c(-Inf, -1, -1)) # <-1 becomes -1

landsatNDVI <- reclassify(landsatNDVI, c(1, Inf, 1)) # >1 becomes 1

3.2 Visualise NDVI

Map out the NDVI values:

plot(landsatNDVI, 
     col = rev(terrain.colors(10)), 
     main = "Landsat 8 NDVI",
     axes = FALSE, box = FALSE)

Plot histogram of NDVI values:

hist(landsatNDVI,
     main = "Distribution of NDVI values", xlab = "NDVI", 
     xlim = c(-1, 1), breaks = 100, yaxt = 'n')
abline(v=0.2, col="red", lwd=2)
abline(v=0, col="red", lwd=2)

3.3 Classify using NDVI

Create a matrix to be used as an argument in the reclassify() function:

reclass_m <- matrix(c(-Inf, 0, 1, #water
                      0, 0.2, 2, #urban
                      0.2, Inf, 3), #veg
                    ncol = 3, byrow = TRUE)
reclass_m
##      [,1] [,2] [,3]
## [1,] -Inf  0.0    1
## [2,]  0.0  0.2    2
## [3,]  0.2  Inf    3

Classify land cover using the defined threshold values:

landsatCover <- reclassify(landsatNDVI, reclass_m)

Plot the land cover classes:


3.4 Save raster

Save the reclassified raster landsatcover in the GeoTiff format:

writeRaster(landsatCover, 
            filename = "clean_data/landsat_landcover.tif", 
            overwrite = TRUE)

4 Landscape metrics

4.1 Quick visualisation

library(landscapemetrics)
library(landscapetools)

#landsatCover <- raster('clean/landsat_landcover.tif') #reload saved raster

show_landscape(landsatCover, discrete = TRUE)

Check if the raster data is in the right format:

check_landscape(landsatCover)
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         3  ✓

4.2 Patch-level

E.g. Area of each patch in the landscape:

lsm_p_area(landsatCover)
## # A tibble: 12,580 x 6
##    layer level class    id metric value
##    <int> <chr> <int> <int> <chr>  <dbl>
##  1     1 patch     1     1 area   14.6 
##  2     1 patch     1     2 area    0.09
##  3     1 patch     1     3 area   27.3 
##  4     1 patch     1     4 area    0.09
##  5     1 patch     1     5 area    0.18
##  6     1 patch     1     6 area    0.36
##  7     1 patch     1     7 area    0.09
##  8     1 patch     1     8 area    0.09
##  9     1 patch     1     9 area    0.09
## 10     1 patch     1    10 area    0.09
## # … with 12,570 more rows

4.3 Class-level

E.g. For each class, the total area of all patches:

lsm_c_ca(landsatCover)
## # A tibble: 3 x 6
##   layer level class    id metric  value
##   <int> <chr> <int> <int> <chr>   <dbl>
## 1     1 class     1    NA ca      6742.
## 2     1 class     2    NA ca     28496.
## 3     1 class     3    NA ca     42933.

E.g. For each class, the average area of patches:

lsm_c_area_mn(landsatCover)
## # A tibble: 3 x 6
##   layer level class    id metric  value
##   <int> <chr> <int> <int> <chr>   <dbl>
## 1     1 class     1    NA area_mn  3.22
## 2     1 class     2    NA area_mn  4.52
## 3     1 class     3    NA area_mn 10.3

4.4 Landscape-level

E.g. Total area of the landscape (all three land cover classes):

lsm_l_ta(landsatCover)
## # A tibble: 1 x 6
##   layer level     class    id metric  value
##   <int> <chr>     <int> <int> <chr>   <dbl>
## 1     1 landscape    NA    NA ta     78170.

5 Credits

Spatial data used in this document:


Creative Commons Licence

Copyright (c) 2020 Song, Xiao Ping